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Abstract 

We study numerically the region above the critical temperature of the four 
dimensional Random Field Ising Model. Using a cluster dynamic we measure 
the connected and disconnected magnetic susceptibility and the connected and 
disconnected overlap susceptibility. We use a bimodal distribution of the field 
with ha = 0.35T for all temperatures and a lattice size L — 16. Through a 
least-square fit we determine the critical temperature at which the two suscep- 
tibilities diverge. We also determine the critical exponents 7 and 7. We find the 
magnetic susceptibility and the overlap susceptibility diverge at two different 
temperatures. This is coherent with the existence of a glassy phase above T c . 
Accordingly with other simulations we find 7 = 27. In this case we have a 
scaling theory with two independent critical exponents. 

1 Introduction 

In the last few years the Random Field Ising Model |T[, or RFIM, has attracted a 
lot of attention. Despite great efforts the critical behaviour of the model is still not 
clear. Both numerical and analytical studies have shown that in three dimensions 
at low temperature and sufficiently small field strength there is a transition from 
a disordered phase to a long range ordered phase. This result was first suggested 
by Imry and Ma 0. They consider the possibility that the model could be split 
in clusters of dimension L. Through a direct comparison between ferromagnetic 
and random field energy they found a value of two for the lower critical dimension 
d[. Subsequent arguments |J, based on perturbative expansion led to the result 
that the critical behaviour of the RFIM should be equivalent to that of the Ising 
model in two dimension fewer. This suggested a dimensional reduction of two so 
that, if this should be taken as a general rule, the lower critical dimension should 
be three instead of two. However, there is a rigorous proof, see Imbrie §, that 
the lower critical dimension is two. It should be shown that, for a certain range of 



1 



temperatures, the mean field equation has more than one solution; this is related to 
the fact that this model has a complex free energy landscape. This is essentially the 
reason why the dimensional reduction failed. An accurate numerical investigation 
of the mean field equation has been done by Guagnelli et al. [13] and successively 
by Lancaster et al. [12]. They have found that the mean field equation has more 
than one solution when the correlation length is still finite. In the spin glass [11] 
mean field theory we have a similar situation. It seems reasonable to use, in this 
case too, replica symmetry breaking theory such as that used by Parisi in that 
context. Mezard et al. ||, [J, using RSB techniques and the SCSA approximation 
|1C], have shown the existence of a region above the critical temperature in which 
there should be a "glassy" phase. In this case we have two different values of the 
critical temperature: one called T c , which is the usual critical temperature of a 
ferromagnetic system, and another , called Tj, so that TJ, > T c , at which we have a 
transition from a paramagnetic phase to a "glassy" phase. In section two we first 
discuss the scaling theory of the model and then we introduce the concept of replica 
susceptibility. In section three we give a brief description of the algorithm used and 
then we show our results and conclusions in the last two section. 



2 Theory 

The RFIM is defined by the Hamiltonian 

H r fi m = —J ^ a i a i — /,hjQ~i- (1) 

<i,j> i 

The variable cr, are Ising like spin and hi are independent random variable with 
mean {hi) av = and variance (hf) av = h R . Typical distribution used is the Gaussian 
or bimodal distribution. 

We first discussed the prediction of the scaling theory. Bray and Moore and 
independently, Fisher || |(| have proposed a scaling theory based on the assumption 
of a second order phase transition with a zero temperature fixed point. At T = 
for the correlation length, as usual, we expect a power law behaviour given by 

£ oc r v . 

In this case t is such that 

f — h3. _ (hR\* 
j j 

where [Hr/ J)* is the value of Hr/ J at the fixed point. The other relevant parameters 
are 

h , J. 

where h is an uniform external field. Because for the RFIM the coupling constant 
is not fixed this yields to a change of the energy scale. In this case we obtain the 
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modified hyperscaling relation 

(2-a) = [d-y)u, 
where y is the critical exponent related to J in particular 

J' = by J. 

Another difference with the Ising model is related to the correlation function. 
The presence of means over the quenched field causes the correlation function to 
have two different types of behaviour. We have a connected and a disconnected 
correlation function. 

T 

G con {r) = < a a r > - < a >< a r > = rd _ 2+v g(r/0- ( 2 ) 

X 

G dis (r) = <a ><a r > = gf(r/Q. (3) 

This defines another critical exponent rj. The (. . .) and the (. . .) denotes, respec- 
tively, the thermal average and the average over the different random field configu- 
ration. The other scaling relations are still valid in this case 

a + 2/? + 7 = 2 , 6 = A//3 ... 

In summary, we have eleven critical exponents and eight scaling relation. There 
seems to be a phase transition with three independent exponents. Schwartz and 
Soffer have demonstrated the inequality 

V < 2r]. (4) 

Numerical simulations O] have suggested that (Q) should be fulfilled like an 
equality. In this case we return to a two independent exponents transition. 

If the transition is a spin glass like transition then the correct order parameter 
of the theory is the overlap 

q = W- (5) 

In this case we are not interested in the magnetisation correlation function but 
in the replica correlation function. More in detailed we can define a magnetic sus- 
ceptibility and an overlap susceptibility 

N[(m 2 ) - (m) 2 ] 
N[(m} 2 ] 

N[W[-W] ( 6 ) 
N[(q) 2 }. 



X{m)dis = 

X{q)con = 

x(q)dis = 
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where N is the volume of the lattice and we have a distinction between the connected 
and the disconnected parts. In this paper by means of Monte Carlo simulation we 
measured the four quantities in (||) in the region slightly above the critical temper- 
ature. In this way we are able to make a comparison between the susceptibility 
related to the magnetisation and that related to the overlap. If the three tempera- 
ture transition scheme proposed in |J] is correct then the overlap and the magnetic 
susceptibility must diverge at two different temperatures. 



3 The algorithm 

The algorithm used to do the simulations is a generalisation of the cluster algorithm 



proposed by Wolff [17] for the Ising Model. According to the limited cluster flipped 
algorithm proposed by Newmann and Barkema [pj| we have realised an algorithm 
capable of flipping more than one spin at a time. The algorithm is capable to forming 
clusters with a limited size. A Monte Carlo step consist of the two following points: 

1. Build a cluster. 



• We choose a random site of the lattice. Then we choose, according to a 
certain distribution probability^], the maximum size, R, of the cluster. 

• We add similarly oriented neighbouring spins. If the spin under consid- 
eration is within the allowed distance then we add it with probability 
1 -exp(-2/3J). 

• We repeat the above step until there are no more spins to add to the 
cluster. 



2. Once the cluster is created we attempt to flip the spins inside it. 

The cluster will be flipped with a probability factor proportional to the random field 
and to the number s of spins which might have been added but which are found just 
outside the radius of the cluster. In detail we have 




exp(-2/3Js) ifm<0 

cxp (—/3m) 
cxp (/3m) 



exp(-/3m) exp( _ 2/3 j s) ifm>0 



where m = J2i^c ^i *- As is explained in [16] this algorithm satisfies the conditions 



of ergodicity and detailed balance for the random field model. A Monte Carlo sweep 
is obtained when we have attempted to flip a number of clusters like the volume of 
the system. 



x As is explained in the appropriate choice of this probability distribution is of fundamental 
importance. In this work we have used a power law distribution i.e. P(R) = 1/R a with a = 2. 
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We believe that for models such as the RFIM, this kind of dynamics is capable of 
strongly reducing the problem related to the dynamic slowing down approaching the 
critical temperature. Using a single spin flip dynamic a new configuration is obtained 
when we try to flip all the spins in the lattice. The probability of flipping a spin 
depends on the local fields through a term proportional to Pfu p (h) = exp(— 2/3/ij). 
It is possible that some spins are aligned to a large local field; in this case such spins 
are almost impossible to be flipped. Sometimes if we flip this "pinned" spin, the 
configuration obtained should be more probable than it was. When one of this spins 
is taken as a part of a cluster the effect of the field of such site should be rounded 
by the other fields in the cluster. In this case the procedure that realises the Monte 
Carlo dynamics is much more complicated than that of the Metropolis algorithm. 
Moreover, this kind of dynamics depends much from the contest. It is therefore 
necessary to spend a lot of time in optimising the algorithm. 



4 Numerical results 

Rieger and Young [14, O] have carried out the most extensive Monte Carlo simula- 



tions in three dimensions in order to test the scaling relation validity. Using finite 
size scaling techniques they calculate all the critical exponents both with a Gaussian 
and with a bimodal probability distribution. 

Making use of the cluster algorithm described in section ^ we carried out Monte 
Carlo simulation in order to search for numerical evidence of the existence of a 
"spin glass" like phase transition in the region above the critical temperature for a 
four dimensional lattice. To this end, at each Monte Carlo sweep, for each disorder 
realisation we measured the average magnetisation (m) its square (m 2 ), the average 
overlap (q) and its square (q 2 ); with m = l/NJ^iLi^i an d q = 1/N J2i=i 
where Oi and r, are two generic spins of two replica of the system. In this way 
we can calculate the four quantities given in (^). The two replica are such that 
they have the same realisation of the disorder. They will approach the equilibrium 
following two different markovian processes, in this case they can be considered 
independent. We performed the calculation using a four dimensional lattice of size 
L = 16 with periodic boundary conditions. We have measured the quantities in (|^) 
for 13 different values of the temperature for each realisation of the disorder. Starting 
from a high temperature region we have cooled the system until ten percent of the 
critical temperature. The hardest region to simulate is the one near the critical 
temperature. For these temperatures the system takes a great deal of time to reach 
equilibrium. Near T c up to 250 Monte Carlo sweeps were needed to balance the 
system. At temperatures T sufficiently greater than T c the system balances very 
fast. In contrast, near T c the time needed to balance the system became too much 
long. For this reason if we took the same amount of Monte Carlo sweeps at each 
temperature we would waste time. We have used a range of temperature varying 
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from 8.33 to 5.97 that correspond to 0.12 < p < 0.1675 where (3 = 1/T. 

Starting from (3 = 0.12 we have cooled the system through the following rules 

P = 0.12 + {k- 1)0.005, (7) 
k = 1...8 

P = 0.1575 + (k- 1)0.0025, (8) 
A: = 9 ... 13. 

In this way we have more points near T c . At each temperature we have used a Monte 
Carlo sweep number (#MCs) given by 

#MCs = 6Ak 2 

k = 1 ... 13 

A third of these have been used to balance the system and the rest to take the 
measurements. For k = 13, i.e. near T c , from the 10816 #MCs, the first 3605 
has been used to reach the equilibrium. This is an order greater than the calculated 
equilibration time. In this way, in the hardest region to simulate the thermal average 
there have been performed over more than 500 uncorrelated measurements. The 
average over the disorder has been performed on 850 samples. We need such a large 
number of samples because the quantities in @ are highly non-self-averaging. For 
this reason the error caused by the disordered average is an order greater than that 
given by the thermal average. The random field has been chosen so that hp; = 0.35T 



for different temperatures. As was pointed out in 14] for greater values of Hr/T 
the system is too difficult to balance and when the ratio is too small the system 
degenerates in the ferromagnetic model. 

The quantity (^), near T c are well fitted by some power law of the tempera- 
ture. Both for the the magnetic susceptibiliy and the overlap susceptibility we are 
sufficiently far from T c so we can neglect finite size effect. For the connected and 
disconnected magnetic susceptibility we used the following power law 



X (M) = C(T-T c )-r (g) 



X (M) dis = d(T-T c y 
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In figures (||) l/x( m )con and 1/ \Jx{ m )dis are plotted against the temperature. 
It is clear from the figure that the two exponents 7 and 7 are respectively slightly 
lower than one and two. The results of the least-square fit are such that 

T c = 5.72 ±0.04 7 = 0.94 ±0.02 7 = 1.91 ±0.08. (10) 

The jacknife technique is used for the errors; the data are correlated since the 
same set of random fields is used at each temperature. According previous simulation 



6 




Figure 1: The full line and the dotted line represent the least-square fit results. The full line is the 



plot of\/f(x) where f(x) = 8.9(T - T c )~ 



id the dotted line is the plot of l//i(a;)' 1/ ' 2 - ) where 



fi(x) = 12.1(T-T C )- 



The intersection with the abscissa gives the critical temperature. 



15] the Schwartz inequality M seems to be valid as an equality. Within the error 
bars we have 27 = 7 that is 

rj = 2rj. 



Considering that < Q 2 > — > 1 when T — > 00 we expect the connected overlap 
susceptibility to have a power behaviour given by 



x(Q\ 



B(T - T b )~" + D 



(11) 



Because of the presence of the random field also the disconnected overlap sus- 
ceptibility has a non vanishing term. We expect a power law behaviour given by 



x(Q)dis = B 1 (T-T b )~ u + D 1 . (12) 

We have reported in figure (|2|) the least-square fit result. As for the results in 
(|l0|) the errors are calculated with the jacknife technique. We find 

T b = 5.88 ±0.04 u = 0.6 ±0.1 U = 0.42 ± 0.05 (13) 

The (^) is valid near the critical temperature. It can be argued that the results 
found for the two critical temperatures, may be an artefact of the power-law be- 
haviour used in (p!^). There could be a temperature drift in the non-singular term 
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as the temperature is not too near from T c . To take care of this effect we can add in 
(O) a linear term in the temperature vanishing near T c . We can use a temperature 
dependence given by, 



x(Q) 



din 



B 1 {T - T*)~" + Di [1 + D Q ^-} 



(14) 



If we fix T* we can perform a four parameters fit. When T* = T c the values 
of x 2 obtained is greater than that obtained in the previous fit. If we set T* = 
we recover the results obtained using (|l|). In this case the temperature can be 
neglected and the results are well represented by the power law given in (12) The 
least-square fit results are reported in tabel [j]. 
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Figure 2: Disconnected overlap susceptibility. The full line is the least-square fit result f(x) 
33(T- 5.88)~ ' 42 ' +41. 



5 Conclusion 

From the data analysis the overlap susceptibility and the magnetic susceptibility 
seems to diverge at two different points. It turns out that > T c . If this is 
the case the three transition scheme, obtained through RSB techniques, should be 
correct. A more extensive Monte Carlo simulations should be done near T c , using 
finite size scaling techniques, in order to confirm the result obtained. It should 
be interesting as well to studing the overlap distribution probability in the region 



8 





( 


13 


T* = T c 


T* = T b 


B x 

7 

D x 
D 

x 2 


33 ±2 
0.42 ±0.02 
41 ±2 

2.14 


13.3 ±0.9 
1.09 ±0.04 

65 ±1 
-0.22 ±0.02 
4.76 


32 ±2 
0.43 ±0.02 

43 ± 2 
-0.03 ±0.02 
1.98 



Table 1: In the first column we have the least-square fit results obtained using the power law given 
in (11). The results in colunm two and three are obtained using the temperature dependence given 
in (14)- In the last raw we have the values of x 2 calculated during the fit. The values in culomn 
one and three are almost equal. When T* — Tt we have Do ~ so that we can assume that fil^ ) is 
a good approximation for all the temperatures used. 



under Another result is given from the comparison between 7 and 7. We find 
7 = 27, according to this result we have found one more scaling relation so that the 
independent critical exponents are two instead of three. In four dimension as well, 
the dimensional reduction is far from giving the correct result, in fact 7 = 0.94 is 
very far from 7/4 which is a prediction of dimensional reduction. 
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